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2.  Summary 


The  Radiative  Transfer  (RT)  approach  is  widely  used  in  applications  involving  scattering 
from  layered  random  media  with  rough  interfaces  (F.T.  Ulaby,  R.K.  Moore  and  AK  Fung, 
Microwave  Remote  Sensing,  Vol.  3,  Artech  House,  1986).  Although  this  approach  involves 
approximations,  in  most  applications  they  are  not  explicitly  stated  or  well  understood.  The 
RT  approach  for  random  media  with  non- scattering  boundaries  has  been  well  studied  (M.I. 
Mishchenko,  Appl.  Optics,  41,  7114-34,  2002).  In  contrast  our  problem  has  scattering 
boundaries  which  are  randomly  rough.  In  order  to  better  understand  the  RT  approach  to 
our  problem  we  adopt  a  statistical  wave  approach  and  then  transition  to  the  RT  equations. 
The  geometry  of  our  problem  consists  of  a  multi-layer  discrete  random  medium  with  rough 
boundaries  which  are  planar  on  the  average.  The  random  medium  in  each  layer  consists  of 
a  homogeneous  background  medium  in  which  discrete  scatterers  are  randomly  distributed. 
The  statistical  characteristics  of  the  random  medium  in  each  layer  are  independent  of  each 
other  and  independent  of  the  statistics  describing  the  rough  interfaces.  The  regions  above 
and  below  the  random  medium  stack  are  homogeneous.  Using  the  Greens  functions  of  the 
problem  without  the  volumetric  fluctuations  we  represent  our  problem  as  a  system  of 
integral  equations.  Employing  the  T-matrix  description  we  first  average  with  respect  to 
volumetric  fluctuations  and  then  use  the  Twersky  approximation  to  obtain  a  system  of 
integral  equations.  We  next  average  with  respect  to  surface  fluctuations,  apply  the  weak 
surface  correlation  approximation  and  arrive  at  a  closed  system  of  integral  equations  for 
the  first  and  second  moments  of  the  electric  fields.  We  use  the  Wigner  transforms  to 
translate  the  coherence  functions  to  radiant  intensities,  which  are  the  fundamental 
quantities  in  the  RT  approach.  On  applying  the  quasi-static  field  approximation  we  hence 
arrive  at  a  system  of  equations  identical  to  those  used  in  the  RT  approach.  From  this  study 
we  learn  that  there  are  more  conditions  involved  in  the  RT  approach  than  widely  believed 
to  be  sufficient. 
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3.  Introduction 


The  radiative  transfer  (RT)  theory  is  widely  used  in  remote  sensing  problems 
[14,28,20,15,1,30].  Often  the  model  of  layered  random  medium  with  rough  interfaces  is 
used.  Multiple  scattering  processes  in  this  structure  are  well  represented  by  the  RT 
equations.  Although  quite  successful  in  numerous  applications  in  various  disciplines,  it  is 
known  that  the  RT  approach  involves  approximations.  Often  people  in  the  remote  sensing 
community  are  not  quite  familiar  with  the  approximations  involved  in  the  RT  approach 
and  hence  there  has  been  inappropriate  use  of  the  RT  approach  in  the  literature.  Since  the 
phenomenological  RT  theory  [7,27]  was  first  developed  for  light  scattering  in  planetary 
atmospheres  the  RT  conditions  prevalent  in  the  atmospheric  context  has  been  popularly 
identified  as  sufficient  conditions  for  employing  the  RT  theory.  However,  we  notice  that 
the  RT  theory  has  been  used  for  a  variety  of  different  problems  [18,24,19,26]  in  various 
applications  with  complex  geometries.  It  is  not  clear  whether  the  classical  conditions 
associated  with  the  RT  theory  are  sufficient  in  all  situations.  In  this  paper  we  will  review 
the  approximations  involved  and  clarify  misconceptions.  In  order  to  better  understand  the 
RT  approach  we  employ  the  more  rigorous  statistical  wave  theory  to  the  problem  and 
hence  make  the  transition  to  the  RT  equations.  In  this  process  we  clarify  and  explain  the 
assumptions  or  approximations  involved  in  the  RT  approach.  By  following  this  procedure 
we  found  that  there  are  more  conditions  embedded  in  the  RT  approach  than  widely 
believed  to  be  sufficient.  For  our  study  we  have  considered  a  multi-layer  random  medium 
composed  of  discrete  scatterers  (see  Figure  1).  By  considering  several  special  cases  of 
this  general  problem  we  show  that  the  number  of  conditions  implied  in  the  RT  approach 
reduces  with  simpler  geometries.  Our  conclusions  are  not  just  for  the  case  discrete 
random  media.  We  show  similar  conditions  apply  for  the  random  continuum  as  well. 

The  report  is  organized  as  follows.  First  we  describe  the  geometry  of  the  problem.  Next 
we  give  the  radiative  transfer  approach  to  the  problem.  The  next  section  is  on  the 
statistical  wave  approach  to  the  problem.  This  occupies  the  major  part  of  the  paper.  It 
deals  with  the  derivation  and  analysis  of  the  first  and  second  moments  of  the  wave 
functions.  A  transition  is  next  made  to  RT  equations.  Next  we  enumerate  and  discuss  the 
conditions  implied  in  the  RT  approach. 
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4.  Description  of  the  Problem 


The  geometry  of  the  problem  is  shown  in  Figure  1.  We  have  an  A-layer  random  medium 
stack  with  rough  interfaces  which,  on  the  average,  are  parallel  planes.  Let  sj  be  the 

permittivity  of  the  background  medium,  and  let  sp  be  the  permittivity  of  the  scatterers  in 

the  y'-th  layer.  The  location,  orientation,  and  shape  of  the  scatterers  are  random  functions 
characterizing  the  fluctuations.  We  assume  the  the  volumetric  fluctuations  in  each  layer 
are  statistically  independent  of  each  other.  Let  N}  be  the  number  of  scatterers,  and  let  p 

be  the  density  of  the  scatterers  (number  of  scatterers  per  unit  volume)  of  the  /-the  layer. 
The  permeability  of  all  the  layers  is  that  of  free  space.  The  randomly  rough  interfaces  are 
given  as  z  =  +  Cj(r__ )  •  Then  C ;  arc  zero-mean  isotropic  stationary  random  processes 

independent  of  volumetric  fluctuations  of  the  problem.  Let  zo  =  0,  and  let  dj  be  the 
thickness  of  the  j- th  layer. 

Let  zo  =  0,  and  let  dj  be  the  thickness  of  the  j- th  layer.  The  media  above  and  below  the 
stack  are  homogeneous  with  parameters  s, ,  ko,  and  eN+1 ,  kN+i,  respectively.  This  system 

is  excited  by  a  monochromatic  electromagnetic  plane  wave  and  we  are  interested  in 
formulating  the  resulting  multiple  scattering  process . 


Figure  1.  Geometry  of  the  Problem 
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5.  Radiative  Transfer  Approach 


Multiple  scattering  in  a  complex  environment  is  well  described  by  the  radiative  transfer 
theory.  This  theory  is  not  only  conceptually  simple  but  also  very  efficient.  The 
fundamental  quantity  here  is  the  specific  intensity  I  which  is  governed  by  the  following 
equation  [7,27,10] 

s-  Vl(r,s)  +  yl(r,s)  =  J  V(s,s')l(r,s')dQ.'  (!) 

One  may  regard  this  equation  as  a  statement  of  conservation  of  energy  density  I  which  is 
a  phase-space  quantity  at  position  r  and  direction  s  .  y  is  the  extinction  matrix  which  is 

a  measure  of  loss  of  energy  due  to  scattering  in  other  directions.  P  is  the  phase  matrix 
representing  increase  in  energy  density  due  to  scattering  from  neighbouring  elements.  Q 
is  the  solid  angle  subtended  by  s .  Given  the  statistical  characteristics  of  the  medium  one 
can  calculate  the  phase  function  using  the  single  scattering  theory  for  elements  that 
constitute  the  random  medium  of  the  layer  [31,11,6,17].  The  extinction  matrix  is  hence 
calculated  using  the  optical  theorem.  The  specific  intensity  in  each  layer  is  governed  by 
an  equation  similar  to  (1).  Since  our  layer  problem  has  translational  invariance  in  azimuth 
the  RT  equation  for  the  m-th  layer  takes  the  following  form, 

COS0iI",(Z’'?)  +  ^,I"'(Z’'?)  =  -L  (2) 


where  the  subscript  m  denotes  that  the  quantity  corresponds  to  that  of  the  m-th  layer  and 
0  is  the  elevation  angle  of  s.  This  set  of  RT  equations  is  complemented  by  a  set  of 
boundary  conditions  which  is  in  turn  based  on  energy  conservation  considerations.  To  be 
more  precise,  we  impose  the  condition  that  the  energy  flux  density  at  each  interface  is 
conserved.  This  leads  to  the  following  boundary  conditions  on  the  m-th  interface 


,  (.*.,$)  =  1  ft  s')  It  (zJ')dCl' + j  X„,„.  (s,s')  C,  (z.s')dCl' 


(3) 


The  boundary  conditions  on  the  (m-1)- th  interface  are  given  as 


K(Zn 


(4) 


where  R  and  X  are  the  local  reflection  and  transmission  Mueller  matrices.  To  be 

mn  mn 

more  specific,  RmB  represents  the  reflection  matrix  of  waves  incident  from  medium  n  on 
the  interface  separating  medium  m  and  medium  n.  The  superscripts  u  and  d  indicate 
whether  the  intensity  corresponds  to  a  wave  travelling  upwards  or  downwards.  These 
Mueller  matrices  are  often  calculated  using  some  asymptotic  theory  such  as  the  Kirchhoff 
approximation  [33,5,30,29].  The  integration  in  these  expressions  are  over  a  solid  angle 
(hemisphere)  corresponding  to  s' .  Suppose  we  have  a  time-harmonic  electromagnetic 
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plane  wave  incident  on  this  stack  from  above.  Then  the  downward  travelling  intensity  in 
Region  0  is 

Io  (z,  s)  =  B0A(cos  0O  -  cos  6,  )S(</>0  -  fa)  (5) 

where  B0  is  the  intensity  of  the  incident  plane  wave  and  {  Oi ,  <7? }  describes  its  direction. 
Since  there  is  no  source  or  scatterer  in  Region  N+l, 

luNJz,s)  =  0 

Notice  again  that  these  boundary  conditions  represent  conservation  of  intensity  at  the 
interfaces.  We  should  point  out  that  the  radiative  transfer  approach  as  applied  to  a 
particular  problem  is  only  a  model  based  on  certain  assumptions.  Since  the  RT  theory  is 
used  in  a  variety  of  applications,  the  particular  assumptions  involved  are  described  in 
terms  of  different  terminologies,  specific  to  the  discipline  where  it  is  used.  One  good  way 
to  understand  in  more  general  terms  the  RT  approach  and  the  underlying  assumptions  is 
to  compare  it  with  the  more  rigorous  wave  approach.  For  the  case  of  an  unbounded 
random  medium  this  kind  of  study  was  carried  out  in  the  1970s  [3,2],  From  that  study  we 
learn  that  the  radiative  transfer  theory  can  be  applied  under  the  following  conditions: 

1 .  Quasi-stationary  field  approximation 

2.  Sparse  distribution 

3.  Statistical  homogeneity  of  the  medium  fluctuations 

These  are  the  well-known  conditions  that  we  associate  with  the  RT  approach.  However, 
our  problem  has  bounded  structures  and,  further,  they  are  randomly  rough.  The  question 
is  this:  are  the  above  conditions  sufficient  to  apply  the  RT  approach  for  our  problem? 
This  is  the  motivation  for  this  paper.  We  follow  the  wave  approach  to  this  problem, 
derive  the  equations  for  the  intensities,  and  hence  make  the  transition  to  the  RT 
equations.  This  procedure  enables  us  to  better  understand  the  necessary  conditions  for 
using  the  RT  approach  for  our  problem. 
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6.  Statistical  Wave  Approach 

The  following  are  the  equations  that  govern  the  waves  in  the  layer  structure: 

Vx  VxE;.  —k'jEj  =  VjEj  j  — 
where 

i= 1 


(6) 

(7a) 


QA  r)  = 


r  e  V,, 


otherwise 


(7b) 


k2  =  m1  he, ,  k]  =  (o2M£j 


(7c) 


JT.  is  the  domain  of  the  7-th  scatterer  and  AT  is  the  number  of  scatterers  in  layer  j.  Thus 
v j  represents  the  volumetric  fluctuations  in  Region  j.  For  the  homogeneous  regions 
above  and  below  we  have 


V  x  V  x  E0  -  E0  =  0 


(8a) 


V  x  V  x  Ew+1  kN+lEN+l  -  0 


(8b) 


The  boundary  conditions  at  the  j- th  interface  are 

n  x  Ej  (r±  ,C/)  =  nx  Ey+1  (r± ,  Cj ) 
n  x  V  x  E ,  (rls  C)  =  n  x  V  x  E  (r±,C,) 


(9a) 


(9b) 


where  n  is  the  unit  vector  normal  to  the  j- th  interface  with  normal  pointing  into  the 
medium  j.  This  system  is  complemented  by  the  radiation  conditions  well  away  from  the 
stack.  We  assume  that  we  know  the  solution  to  the  problem  without  volumetric 
fluctuations,  and  denote  it  as  E  .  Let  the  Green's  functions  to  this  problem  be  denoted  as 


These  Green's  functions  are  governed  by  the  following  set  of  equations: 
V  x  V  x  6 (r,  r 1 ’)  -  kfi  ,  (r,  r')  =  lSlkS(r  -  r ') 


(10a) 
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(10b) 


nxGjk  (t±  ,  Cj  ;r')  =  nx  G{j+V)k  (r± ,  £ j ;  r ') 
n x V  x  G  jk  (r±,  £ , ;  r')  =  n  x  V  x  G(  /+1),  (r± ,  £j ;  r') 


(10c) 


where  I  is  unit  dyad.  The  boundary  conditions  above  are  for  the y-th  interface.  Similarly 
on  the  (/-  l)-th  interface  we  have  the  following  boundary  conditions. 


h  x  G  jk  (rx ,  £  ;r')  =  nxG(  H)k  (r± ,  ;  r') 


(11a) 


nxVxG, (r± , ;r')  =  nxVxG,  l)k (r± , ; r') 


(lib) 


Using  these  Green's  functions  and  the  radiation  conditions  the  wave  functions  can  be 
represented  as 


N 

E.(r)  =  E.(r)  +  Xf, 

k= 1 


G,,(r,rX(r')E,(r')Jr' 


(12) 


where  QA  =  {rj;^  <  z'  <  Xi}  •  Note  that  vo  =  vn+ i  =  0  • 


In  order  to  carry  out  multiple  scattering  analysis  with  a  distribution  of  discrete  scatterers 
it  is  convenient  to  employ  the  concept  of  transition  operator  T  [13,34,32].  Suppose  we 
know  the  electric  field  E;  incident  on  the  /-th  scatterer.  We  introduce  the  transition 
operator  T'  such  that  the  electric  field  scattered  by  the  /-th  scatterer  is  given  as  T'E' . 
Using  this  concept  (12)  may  be  expressed  as 


N  Nt 


E,(r)  =  E,(r)  +  m.  dr'ln  dr"G jk(r ,r')Tlk(r' ,r")Elk(r")  j  -  0,1,...,/V  +  1  (13) 


it=l  l=\ 


Note  that  TA  depends  only  on  the  /-th  scatterer.  To  proceed  further  with  the  multiple 
scattering  analysis  it  is  expedient  to  use  symbolic  representation  of  (13). 

E,=VW  (14) 


First  we  average  (14)  w.r.t.  volumetric  fluctuations  to  get 

e/)„=T+G(tX), 


(15) 
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where  the  subscript  v  denotes  volumetric  averaging.  Since  there  are  Nk  scatterers  in  the 
k- th  layer 


tX 


=  J  TX/X ,  r2 , . . . ,  ;  Sj ,  s2 , . . . ,  )  dr  ds 


(16) 


where  p  is  the  joint  probability  density  function  of  finding  the  scatterers  at  r, ,  r2 , . . . ,  r;V, 
with  orientations  s1,s2,...,s2V  We  assume  that  the  positions  and  orientations  are 
independent  of  each  other.  In  other  words 


p{rx,r2,...,r ^  ;s1,s2,...,s„  )  =  />(r„r2,...,r N  )p{s„s2,...,sN  ) 


(17) 


Furthermore,  assume  that  the  orientation  of  the  particle  at  position  r,  is  independent  of 
the  orientation  of  all  other  particles,  which  means 


p(sns2,. 


sX  =  f>(s,) 


j= 1 


(18) 


We  next  express  the  joint  position  probability  density  function  as 
P(rl,r2,...,rlfk)  =  p(rl,r2,...,r'l,...,rNk  ^^(r,) 


(19) 


where  p(rt ,  r2 , ....  r/ , ... ,  rVj  |  r, )  is  the  conditional  pdf  of  finding  scatterers  at  r, ,  r2 , . . . ,  rv. 

by  fixing  the  /-th  scatterer  is  at  r, .  The  prime  in  r/  denotes  that  r,  should  be  omitted  in 

the  argument  list  of  the  conditional  probability  density  function.  Substituting  this  relation 
in  (16)  we  obtain 


(20) 


where  Ej  'j  denotes  conditional  average  with  scatterer  /  fixed  at  r, ,  If  W  is  large  and  the 
distance  between  scatterers  is  large  then  we  can  approximate 

(e;  )'□(£;)  <2» 

This  is  called  the  Foldy’s  approximation  and  is  applicable  for  sparse  media.  Under  this 
approximation  (15)  becomes 

(E/Uf  +/>A.(T.)<E-).  <22) 
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where  pk  is  the  density  of  the  scatterers  in  layer  k.  Now  operate  the  above  equation  by 
VxVxI-^2I  to  obtain 


VxVx(E,.)i-t;(E() 


V 


(23) 


Next  average  this  over  surface  fluctuations  to  get 


VxVx(E;)  -k 


T,  E 


(24) 


Note  that  the  transition  operators  are  independent  of  surface  fluctuations.  From  this  we 
see  that 


VxVx/eA  -kU E.\  =0  j  =  0,N  +  l 

\  J  l  VS  J  \  J  /  vs 


(25) 


which  means  that  the  coherent  propagation  constants  in  the  regions  above  and  below  the 
layered  stack  are  unaffected  by  the  fluctuations  of  the  problem.  However  they  indeed  get 
modified  in  the  layered  stack  region.  On  writing  (24)  as 


{VxVxI-*;-/>,(Ti}}(E,}jj=0 


(26) 


we  infer  that  x,  =  yjk/  +  P,  (t,-)  represents  the  mean  propagation  constant,  in  operator 
form,  for  coherent  waves  in  layer  j. 

Since  the  problem  is  statistically  homogeneous  in  azimuth  the  mean  fields  in  our  system 
have  the  following  form: 


E /(i0)  =exp(iku-r) 


A-  (kXi)p+  exp 


iq^z 


+  5f(k±;)p  exp 


-WjZ 


7=1,2, 


,N 


(27) 


(Eo  (r)}vs  =  exp(/k±;  •  r)  {p'  exp  [-ik^z]  +  Rp  (k1;.  )p+0  exp  [z^.z]}  ^ 

{K+l  (r))vs  =  exp(/k±i  •  r)Xp  (k1()p,,+l  exp  \_-ikp(N+l)ziz\  (2 

where  the  superscript  p  stands  for  the  polarization,  either  horizontal  or  vertical,  p  is  the 
unit  vector  representing  polarization.  <7  4s  the  z-component  of  /  ,  The  subscript  i  is  used 

to  indicate  that  the  wave  vector  is  in  the  incident  direction.  R  and  X  denote  respectively 
the  mean  reflection  and  transmission  coefficient  of  the  stack.  and  B .  denote 

respectively  the  mean  coefficients  of  up-going  and  down-going  waves  in  the y-th  layer. 
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Based  on  this  we  can  formulate  the  waves  averaged  w.r.t.  volumetric  fluctuations  as 


E;(0),~j  expO'k, .r){y«(ki,kil>;  e"v ,  (30) 

Eo(r))y  =  exp(/k  ,  •  r)exp[-/7c0r;z]p0  +~^rj  exp(/k  ■r)R''q(kl  M  i)q;  exp [ik0zz]dU  (3]) 


and 


Ev+1(r))v  =^r{ exp(ik± -r)Xpq(k±,ku)q-N+l  exp \_-ik(N+l)zz\  dk 


(32) 


where  A.  ,Bf,R  and  X  are  now  integral  operators  representing  scattering  from  rough 

interfaces.  The  boundary  conditions  associated  with  the  above  equations  at  the y-th 
interface  are 

hx(E.(r1,^)}i,=nx(E.+1(r±,^))i, 


j  =  1,2, ...,7V 


(33) 


and 

n<Vx(EJ(r1,fJ))i  =  iixVx{EJ,,(r1,f,})v  j=l,2,...,N 

The  above  system  may  be  solved  either  numerically  or  by  any  one  of  the  asymptotic 
methods  available  in  rough  surface  scattering  theory  [5,4,33]  to  evaluate  the  mean 
coefficients  that  appear  in  (27)-(29). 


We  proceed  now  to  the  analysis  of  the  second  moments,  by  starting  with  (12).  For 
convenience  we  write  it  in  symbolic  form  as 


E;=E, 


+  ZC#T,'E 

k= 1 


(35) 


We  take  the  tensor  product  of  this  equation  with  its  complex  conjugate  and  average  w.r.t. 
volumetric  fluctuations  and  obtain 


Ey  ®E/)  =  (E/ 


E 


N  N  Nk  Nk. 

+  IIIZ  (G* 

k= 1  k'= 1  1=1  l'= 1 


Gy  KwE;®Et 


(36) 
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where  K  is  the  intensity  operator  of  the  volumetric  fluctuations.  Employing  the  weak 
fluctuation  approximation  we  approximate  K  by  its  leading  term 


K„n(T'le>Tll-)sa.s„.\ 


(37) 


On  substituting  this  in  (36)  we  get, 


E- ®E‘>  =(E,>  ®(E-) +Z  ft (Gft)  ®(G;)„(Tt ®t;)/e1  ®e; 


k= 1 


Next  we  average  (38)  w.r.t.  the  surface  fluctuations  to  get 


j  j 


N 


=((E()  ®(EU  )  +X  Pk  ((Gft )®  c;,)  }  (Tk®Tk)(Ek®Ek 


k= 1 


where  we  have  used  the  following  approximation 

A), . ®(g;),(t» ®t;)(e,  ®e;)  )  □  ({G„\ (T,  ®t;){e,  ®e; 


(38) 


(39) 


(40) 


We  call  this  the  ‘weak  surface  correlation’  approximation,  which  we  will  see  later  to  be 
an  important  condition  embedded  in  the  RT  approach  to  our  problem. 

As  it  stands  this  equation  is  very  difficult  to  solve  either  analytically  or  numerically. 
Besides,  one  important  goal  for  us  is  to  investigate  the  conditions  needed  for  employing 
the  radiative  transfer  approach  for  our  problem.  With  these  in  mind  we  introduce  Wigner 
transforms.  Note  that  (39)  is  an  equation  for  the  coherence  function  which  is  a  ‘space- 
space’  quantity.  On  the  other  hand  the  RT  equation,  as  we  saw  earlier,  is  an  equation  for 
the  specific  intensity  which  is  a  ‘phase-space’  quantity.  Wigner  transforms  serve  as  a 
bridge  to  link  these  two  quantities  [35,8,16,22], 


We  introduce  Wigner  transforms  of  waves  and  Green's  functions  as 


r +  r 


'  k  =  J  (E™  (r  ))v  0  (E'«  (r'))vs  eXP  [~ik  •  (r  -  r')]  ^(r  -  r') 


(41a) 


fr  +  r' 

l  2 


’ k  r  1  ((E™  (r  ))v  0  (E™  (r')))s  exp  l~ik  •  (r  - r')]  d(r  - r') 


(41b) 


(  r  +  r'  .  ,  r,  +  r,' 


k  I  -  V.l  H  d{Y  -r,)l  J(r>  ~0^r,)  e'|■(r,_r,',  \{^rnn  (r  ’  r!  ))v  ®  (Glm  (r ',  r/ 


(42) 
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(43) 


In  terms  of  these  transfonns  (39)  becomes 


A  1 V 

^m(r’k)  =  ^m(r’k)  +  ^7rX^»  J  dr'^da^d$3mn(r,k\r',a)Tn(r',u\r",p)£n{r",p) 


n= 1 


with 


T„(R15a|  R2,P)  =  Jjr,Jjr2  exp  {-/a  ■  r,  +  /p  •  r2 }  /t(r,  +  -|,R2  +  -!|-)®T*(ri  -y,R,  -y) 


(44) 


where  T  is  the  element  transition  operator  in  n-th  layer. 


The  fact  that  our  problem  has  translational  invariance  in  azimuth  implies  the  following: 


(r5k)  =  6m  (z,k)  (45a) 

(r>k  I  r',l)  =  0mn  (z,k  I  z',  1; r±  —  r[ )  (45b) 

Using  these  relations  in  (43)  we  have 

Sm  (Z’  k)  =  K  iZ’k)  +  tAtZ  Pn  Jr  dz'\  da\ d^mn  (Z’ R  I  Z'’  «’  0K  (<*>  P)S »  i2' ’  P) 

\Z7l )  n=\  » 

where 

smn  (z,  k  I  z\  a;  0)  =  J  3mn  (z,  k  |  z',  a;  r±  -  r{  >/(r±  -  r[ ) 


(46) 

(47) 


eT(a,/?)  =  f(a,/?)0f‘(a,/?)  (48) 

f  is  the  element  scattering  matrix  in  the  n-th  layer.  Since  the  medium  is  assumed  to  be 
sparse  inter-particle  scattering  takes  place  in  the  far-field  zone  of  each  other.  It  is  based 
on  this  fact  that  we  have  transitioned  from  Tn  to  JF  .  Also  we  have  employed  the  on-shell 

approximation  to  Tn . 

To  proceed  further  we  need  to  evaluate^.  Further  we  need  to  relate  this  system  with 
that  of  RT,  which  involves  the  boundary  conditions  at  the  interfaces.  Therefore  we  need 
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to  identify  the  coherence  functions  corresponding  to  up-  and  down-going  wave  functions. 

To  facilitate  this  we  decompose  (Gmn)  into  its  components, 

(G  )  =5  G°  +G““  +G“rf  +Gdu  +  Gdd  (49) 

\  mn  I  v  mn  m  mn  mn  mn  mn 


where  the  first  tenn  is  the  singular  part  of  the  Green's  function.  The  superscripts  u  and  d 
indicate  up-  and  down-going  elements  of  the  waves.  The  other  components  are  due  to 
reflections  from  boundaries.  These  are  formally  constructed  using  the  concept  of  surface 
scattering  operators  as  follows  [33], 


G*  (r,r'f  = 


(2^T 


J <fk±J d\t'L  {S'“*(k1,k'L)},m  exp |/k±  -r +  m^'(k1)z}exp{-ik;  ■r,-ibqvn(V!jz'}  (50) 


where  is  the  surface  scattering  operator.  The  superscripts  a  and  b  on  S  are  used  to 
indicate  whether  the  waves  are  up-going  or  down-going.  In  the  exponents,  a,b= 1  if  the 
waves  are  up-going.  We  let  a,  b=- 1  if  the  waves  are  down-going.  The  z-component  of  the 
mean  propagation  constants  in  the  n- th  layer  is  denoted  as  qn  .  We  recall  that  ebnm  is  the 

Wigner  transform  of  ((Gmn)  ® (G^n^  ^  .  The  superscripts  //. v  stand  for  polarization, 

either  h  or  v.  When  we  use  (34)  to  perform  the  Wigner  transfonn  we  ignore  all  cross 
terms.  In  other  words,  we  make  the  following  approximation, 

§  = 8  S°  +SUU  +Sud  +Sdu  +Sdd  (51) 

mn  mn  m  mn  mn  mn  mn 

where  3“bn  is  the  Wigner  transform  of  ®  ^  .  This  approximation  reminds 

us  of  the  concept  of  incoherent  addition  of  intensities  associated  with  the 
phenomenological  radiative  transfer  theory. 


With  the  introduction  of  this  representation  for  Smn  in  (43)  we  can  trace  up-  and  down¬ 
going  waves  to  obtain  the  following  equations  for  the  coherence  function: 


(2 


Pm\[  dz'jdaj (z.k  |  z’,o;0)T”(o, p)S„ 


(2 t 


dz'j  daj (z,k  |  z',a;0)T*(a 


(52a) 


^(z,k)  =  £:l(z,k) 
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(2*y 

4  N 


•/>,};'  dz'j daj JPC (z,k  I  z', a; 0)Y* (a, p)&“ (z', /?) 


(2^)4 


Zp„r  Jz'l  M  <**£  (z’k  I  ^«;0)Tf  {a,P)&bn{z\  p) 


(52b) 


Note  that  summation  over  a,b  =  { u,d }  is  implied  in  the  above  equations.  The  first  term 


in  these  equations,  ssa ,  represents  the  contribution  due  exclusively  to  surface  scattering, 
and  has  the  following  form: 


{E:,r>:;;r'(kJ 


'Li 


(53) 


where  2“  is  the  amplitude  of  the  up-going  wave  in  the  m-th  layer  after  volumetric 
averaging  is  performed.  This  means  that  it  is  a  random  function  of  surface  fluctuations. 
When  we  substitute  (53)  and  the  expressions  for  ~smn  in  (52)  we  find  that 


K(z,k)}^=2xslkz  -|a[^(k±)  +  ^*(k_L)]|exp 


m 


(54) 


On  substituting  this  in  (52)  and  differentiating  w.r.t.  z  we  obtain  the  following  transport 
equations: 


-  *  \ju  (k± )  -  L  (k±  )]|  S;v  (z,  k±)  = 

=  4  P„,  J  da  SI  (a± )  ®  SI*  (a± ) 


jk± 4  (k±)  +  &  (k± )] ;  1 1 a  [$V  («i )  +  £  («±)]1  K'v'  (z’‘ a± ) 


(55a) 


\_d__. 

,  dz 


?„(ki)-sC(k  ±)JJ*v(*>k±)  = 

=  V„J^a±^(ax)®C(«i) 

rkj.’-^[^(ki)  +  ^(k±)];«i4fl[^'(a^)  +  ^'(a^)]V/'v(^«i) 


, , 

Hv\li  v 


(55b) 


where  summation  over  a  is  implied.  When  the  superscript  a  corresponds  to  u  the  value  of 
a  in  the  argument  of  ym  take  the  value  +1;  on  the  other  hand  when  the  superscript  a 
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corresponds  to  d  the  value  of  a  in  the  argument  of  ym  take  the  value  -1.  Since  all 
quantities  in  (55)  correspond  to  the  same  layer  m  we  have  dropped  the  subscript  in  in  y 
and  s  to  avoid  cumbersome  notations.  To  obtain  appropriate  boundary  conditions  we 
have  to  go  back  to  the  integral  equation  representations  for  eu  and  sd ,  examine  their 

behavior  at  the  interfaces,  and  seek  a  relation  between  them.  After  considerable  effort  we 
managed  to  arrive  at  the  following  boundary  conditions.  At  the  (m-  l)-th  interface  we 
have 


(56) 


with  tH  =  R®R  where  Rm_lm  is  the  stack  reflection  matrix  (not  the  local  reflection 
matrix)  for  a  wave  incident  from  below  on  the  (m-1)- th  interface.  Similarly 

C(z„-„k1)  =  |(3i„,.„(k±,k;)}sr,J(z„.„k;)rfk; 


(57) 


where  9?m+1  m  is  the  tensor  product  of  stack  reflection  matrix  for  a  wave  incident  from 

above  on  the  (m-l)-th  interface.  We  were  able  to  obtain  the  boundary  conditions  only 
after  imposing  certain  approximations  such  as  the  one  given  below.  Consider  the 
following  identity: 


cdu  =  n  r  (c>  +suu  1 T) 

^ mm  /«  vm— l,m  m  ^ mm  j  n 


(58) 


where  Dm  =  di ag  { exp(iqhdm ),  exp(iqvdm  )}  .  Notice  that  this  is  an  operator  relation  where 

all  elements  are  operators.  Taking  the  tensor  product  of  (58)  with  its  complex  conjugate 
we  have 


-i  du 
^  mm 


Sd“*  =  (D  ®D*)  R 

mm  V  m  m  '  \ 


(x)  R 


{s 


,  (Tt  UU 

m  ^  mm 


){S 


,  r\llU 

m 


(D  ®D*  ) 

\  m  m  J 


(59) 


Next  we  average  (59)  w.r.t.  surface  fluctuations  and  get 


s';  ®  sr)  □  (d.  ®  d;  )  (r„_„  ®  r;.,, 

{s;+si}®{s;+s; 


(D  ®D*  ) 

\  m  m  J 


(60) 


where  we  have  approximated  that  the  two  tensor  products  in  the  middle  are  weakly 
correlated.  A  further  approximation  that  we  impose  is  given  as  follows: 


(s>  +s"u  }®{s>  +s™  Y\u  s>  ®s>*+(s"u 

[  m  mm  J  m  mm  J  /  m  m  \  mm  mm  j 


(61) 


These  are  the  kinds  of  approximations  required  to  arrive  at  our  boundary  conditions. 
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7.  Transition  to  Radiative  Transfer 


Next  we  have  to  transition  from  this  transport  equation  (55)  to  the  phenomenological 
radiative  transfer  equation  discussed  earlier.  To  accomplish  this  we  have  to  link  the  key 
quantities  of  waves  and  radiative  transfer,  viz.,  coherence  function  and  specific  intensity. 
The  relation  between  them  is  obtained  by  computing  the  energy  density  using  the  two 
concepts.  Thus  we  have 

|  e  (Em  (r  )E\  (r))  =  JJw  J  Iuv  (r ,  s)d£\ 


(62) 


Wigner  transform  provides  us  following  relation: 

(£„(r)£») 


(63) 


From  (62)  and  (63)  we  relate  /  to  s  as 

1  k'2  (64) 

/-<z-S)=^^cos&-(z^> 

Now  we  can  transition  to  the  phenomenological  RT  equations.  Using  the  relation 
between  &  and  /  we  change  the  integration  variable  to  solid  angle  and  arrive  at  the 
following  equation, 

|cos^+ n. }  /;  (*,  •?)  =  { [p;u  (n,  O')  Ij  (z,  s) + p;“  (n,  n  ')  ^  (z,  S')}dn' 

|- cos  e  j- + Vij /;  (z,  s)  =  J  {if  (q,  Q' ■)/;  (z„ s') + P«d  (Q,! Q')  Idj  (z„ 

where  77  is  the  extinction  matrix  and  P  is  the  phase  matrix.  Implicit  summation  over 
subscript  j  is  assumed  in  (65).  To  facilitate  comparison  with  the  results  of  Ulaby  et  al. 

[30],  and  Lam  and  Ishimaru  [12]  we  have  used  a  modified  version  of  Stokes  vector  [10]. 

Instead  of  the  standard  form  {/,  Q,  U,  V}  we  use  {(I  +  Q)/2,(I -Q)/2,U,V}  .  The 

subscript  of  /  denotes  the  element  number  of  our  modified  Stokes  vector.  Although  the 
structure  of  this  equation  is  identical  to  that  of  the  RT  (equation  (2)),  the  elements  of  the 
phase  matrix  and  the  extinction  matrices  are  not  the  same.  The  primary  reason  is  because 
of  the  differences  in  the  real  part  of  the  mean  propagation  constants  of  horizontally  and 
vertically  polarized  waves.  On  assuming  that  q'h  =  q'v  =  k'mz  we  obtain  the  following 
expressions  for  the  extinction  and  phase  matrices: 

=  -  cos  6>  diag  {2q"v,  2 q"h ,  q"v  +  q"h ,  q"v  +  q"h }  (66) 


(65a) 

(65b) 
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7f(Q,Q')  =  £?.- 

( k± ,  ak  cos  0;  k^ ,  M  cos  O' } 

(67) 

where 

II 

_  ri 

%  =(Re(/„A)> 

i?4=-( 

tm  (/„./!)) 

II 

< 

K> 

% =(|/J2} 

%  =(Re(/„/;)) 
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We  have  suppressed  the  arguments  for  brevity.  For  instance, 
{k1,akcos0;k'L,hkcos0'} 

=  (Re  {  f vv(k±,ak  cos  d;k'±,bk  cos  d')f*h  (k±,<2kcos0;k'L,6kcost?')}^ 


Note  that  these  transport  equations  (65)  are  identical  to  those  of  classical  RT  equations 
(2)  that  we  described  in  Section  2.  Thanks  to  our  statistical  wave  approach  we  now  have 
explicit  expressions  for  the  extinction  matrix  and  phase  matrix  in  terms  of  the  statistical 
parameters  of  the  problem.  Let  us  now  next  turn  our  attention  to  the  boundary  conditions 

(BC).  In  our  wave  approach  we  obtained  BCs  in  terms  of  'stack'  reflection  matrix  R, 
whereas  in  the  RT  approach  the  BCs  are  given  in  terms  of  the  local  interface  reflection 
matrices.  We  can  readily  reconcile  with  this  apparent  difference.  Note  that  the  BC  in  the 
wave  approach  forms  a  closed  system  whereas  in  the  RT  approach  it  is  open  (linked  to 

adjacent  layer  intensities).  Let  us  take  a  look  at  the  BC  at  the  (m-l)-th  interface.  Rm_2m 

can  be  expressed  in  terms  of  Rm_2  m_,  as  follows, 
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! — 1,1 
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-i,/ 


+  T,„ 


-i-i 


I-R 


,D„ 
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i-2,i 


,D  _jT 


(70) 


This  is  the  relation  between  the  stack  reflection  coefficients  of  adjacent  interfaces.  The  R 
and  T  are  local  (single  interface)  reflection  and  transmission  matrices  at  the  (m-l)-th 
interface.  On  operating  E"  with  (70)  we  get 
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(71) 


E 


d 

m 


R  .  E  +T  ,E 

m—l,m  m  m,m—\ 


d 

m- 1 


Notice  that  this  boundary  condition  now  involves  only  local  interface  Fresnel 
coefficients.  Take  the  tensor  product  of  (71)  with  its  complex  conjugate  and  average 
w.r.t.  surface  fluctuations.  Employing  the  Wigner  transform  operator  on  this,  we  obtain  a 
boundary  condition  at  the  (m- l)-th  interface  similar  to  that  of  the  RT  system.  However, 
the  reflection  and  transmission  matrices  used  in  the  RT  system  correspond  to  unperturbed 
medium  as  opposed  to  the  average  medium  as  in  the  case  of  the  wave  approach. 


Similarly  we  write  Rm+lm  in  terms  of  Rm+2m+1  and  hence  obtain  the  BC  at  the  m-th 


interface  as 


E"  =  R  Erf  +  T  E" 

m  m+\,m  m  m,m+ 1  m+ 1 


(72) 


Take  the  tensor  product  of  (71)  with  its  complex  conjugate  and  average  w.r.t.  surface 
fluctuations.  Employing  the  Wigner  transfonn  operator  on  this,  we  obtain  boundary 
condition  at  the  m-th  interface  identical  to  (3)  (after  making  the  approximation  as  before). 
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8.  Discussion 


Now  that  we  have  made  the  transition  from  statistical  wave  theory  to  radiative  transfer 
theory  it  is  instructive  to  itemize  the  assumptions  implicitly  involved  in  the  RT  approach: 

1 .  Quasi-stationary  field  approximation 

2.  Sparse  distribution  of  scatterers 

These  are  the  well-known  conditions  necessary  for  the  unbounded  random  medium 
problem.  However,  if  the  medium  is  bounded  we  need  to  impose  additional  conditions. 
We  found  that  the  extinction  coefficients  calculated  in  the  wave  approach  and  the  RT 
approach  are  different  and  only  after  applying  further  approximations  can  they  be  made 
to  agree  with  each  other.  The  following  additional  condition  is  required  for  our  bounded 
random  medium  problem: 

3.  Layer  thickness  must  be  of  the  same  order  or  greater  than  the  corresponding  mean 
free  path. 

When  the  interfaces  are  randomly  rough  we  need  the  following  additional  conditions. 

4.  Weak  surface  correlation  approximation. 

5.  All  fluctuations  of  the  problem  are  statistically  independent  of  each  other. 

Recently  Mishchenko  et  al.  [17]  (hereafter  referred  to  as  MTL  for  brevity)  derived  the 
vector  radiative  transfer  equation  (VRTE)  for  a  bounded  discrete  random  medium  using  a 
rigorous  microphysical  approach.  This  enabled  them  to  identify  the  following 
assumptions  embedded  in  the  VRTE. 

1 .  Scattering  medium  is  illuminated  by  a  plane  wave. 

2.  Each  particle  is  located  in  the  far  field  zone  of  all  other  particles  and  the 
observation  point  is  also  located  in  the  far  field  zones  of  all  the  particles  forming 
the  scattering  medium. 

3.  Neglect  all  scattering  paths  going  through  a  particle  two  or  more  times  (Twersky 
approximation). 

4.  Assume  that  the  scattering  system  is  ergodic  and  averaging  over  time  can  be 
replaced  by  averaging  over  particle  positions  and  states. 

5.  Assume  that  (i)  the  position  and  state  of  each  particle  are  statistically  independent 
of  each  other  and  those  of  all  other  particles  and  (ii)  spatial  distribution  of  the 
particles  throughout  the  medium  is  random  and  statistically  uniform. 

6.  Assume  that  the  scattering  medium  is  convex. 

7.  Assume  that  the  number  of  particles  N  forming  the  scattering  medium  is  very 
large. 

8.  Ignore  all  the  diagrams  with  crossing  connections  in  the  diagrammatic  expansion 
of  the  coherency  dyadic. 

Below  we  make  a  connection  between  the  two  by  considering  each  condition  derived  by 
MTL  and  relating  it  to  ours.  We  will  denote  the  condition  numbers  derived  by  MTL  as 
MTL  #  and  those  obtained  in  this  paper  as  SM  # 

•  MTL  1:-  We  also  have  a  plane  electromagnetic  wave  illuminating  our  system, 
although  as  pointed  out  by  MTL  it  can  be  a  quasi-plane  wave. 
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•  MTL  2:-  We  also  have  implicitly  employed  the  far  field  approximation.  It  is 
embedded  in  SM  1  and  SM  2. 

•  MTL  3:-  This  is  embedded  in  SM  2.  Although  not  explicitly  stated,  the  scattering 
processes  as  mentioned  in  MTL  3  are  avoided.  This  is  identical  to  the  Foldy's 
approximation  use  in  this  paper. 

•  MTL  4:-  In  our  paper  we  have  restricted  our  attention  to  the  time-independent 
problem  and  hence  did  not  encounter  the  issue  of  ergodicity. 

•  MTL  5:-  We  have  also  employed  this  assumption  although  we  did  not  itemize 
this  as  a  condition. 

•  MTL  6:-  In  our  problem  we  have  distinct  scattering  boundaries  and  the  character 
of  the  waves  exiting  or  entering  them  are  explicitly  contained  in  the  boundary 
conditions.  Hence  convexity  of  the  scattering  medium  is  not  a  necessary  condition 
for  us. 

•  MTL  7:-  This  condition  is  embedded  in  SM  3. 

•  MTL  8:-  This  condition  is  embedded  in  SM  2.  Under  sparse  medium  conditions 
we  only  take  into  consideration  the  leading  term  of  the  intensity  operator. 

Since  the  problem  we  considered  in  this  report  involves  scattering  boundaries  we  have 
some  additional  conditions  beyond  those  of  MTL.  There  are  a  few  more  remarks  that  we 
would  like  to  make  before  closing,  (a)  In  RT  theory  the  medium  is  assumed  to  be  sparse 
and  hence  the  “refraction  effects”  of  the  fluctuations  are  ignored.  Thus  in  the  boundary 
conditions  we  should  use  the  background  medium  parameters  rather  than  the  effective 
medium  parameters  as  derived  in  our  statistical  wave  theory,  (b)  To  arrive  at  (50)  from 
(49)  we  have  ignored  the  contribution  of  evanescent  modes. 

To  summarize,  we  have  enquired  into  the  assumptions  involved  in  adopting  the  radiative 
transfer  approach  to  scattering  from  layered  random  media  with  rough  interfaces.  To 
facilitate  this  enquiry  we  adopted  a  wave  approach  to  this  problem  and  derived  the 
governing  equations  for  the  first  and  second  moments  of  the  wave  fields.  We  employed 
Wigner  transforms  and  transitioned  to  the  system  corresponding  to  that  of  radiative 
transfer  approach.  In  this  process  we  found  that  there  are  more  conditions  implicitly 
involved  in  the  RT  approach  to  this  problem  than  it  is  widely  believed  to  be  sufficient. 
With  the  recent  development  of  fast  and  efficient  algorithms  for  scattering  computations 
and  the  enormous  increase  in  computer  resources  it  is  now  feasible  to  take  an  entirely 
numerical  approach  to  this  problem  without  imposing  any  approximations. 

In  spite  of  such  developments,  to  keep  the  size  of  the  problem  manageable  only  special 
cases  have  been  studied  thus  far  [9,23,21,25].  Hence  it  is  very  much  of  relevance,  interest 
and  convenience  to  apply  the  RT  approach  to  these  problems.  However,  one  should  keep 
in  mind  the  assumptions  involved  in  such  an  approach.  Otherwise  interpretations  of 
results  based  on  RT  theory  can  be  misleading. 
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